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ABSTRACT 

By performing a global magneto-hydrodynamical simulation for the Milky Way with 
an axisymmetric gravitational potential, we propose that spatially dependent ampli¬ 
fication of magnetic fields possibly explains the observed noncircular motion of the 
gas in the Galactic center region. The radial distribution of the rotation frequency 
in the bulge region is not monotonic in general. The amplification of the magnetic 
field is enhanced in regions with stronger differential rotation, because magnetoro- 
tational instability and field-line stretching are more effective. The strength of the 
amplified magnetic field reaches 0.5 mG, and radial flows of the gas are excited by 
the inhomogeneous transport of angular momentum through turbulent magnetic field 
that is amplified in a spatially dependent manner. In addition, the magnetic pressure- 
gradient force also drives radial flows in a similar manner. As a result, the simulated 
position-velocity diagram exhibits a time-dependent asymmetric parallelogram-shape 
owing to the intermittency of the magnetic turbulence; the present model provides 
a viable alternative to the bar-potential-driven model for the parallelogram-shape of 
the central molecular zone. This is a natural extension into the central few 100 pc of 
the magnetic activity, which is observed as molecular loops at radii from a few 100 pc 
to 1 kpc. Furthermore, the time-averaged net gas flow is directed outward, whereas 
the flows are highly time-dependent, which we discuss from a viewpoint of the outflow 
from the bulge. 

Key words: accretion, accretion disks — Galaxy: bulge — Galaxy: centre — Galaxy: 
kinematics and dynamics — magnetohydrodynamics (MHD) — turbulence 


1 INTRODUCTION 

Observations of the atomic and molecular interstellar 
medium in the inner Milky Way exhibit large noncircular 


& Burton] 19781 

Bally et al. 

1987 

Oka et al.||1998| |Tsuboi, 

Handa & Ukita| 

1999| ISawac 

la et al.||20041 Oka et al.||2005| 


for the noncircular motion, de Vaucouleurs (1964) raised a 


possibility that the gas responds to a stellar bar potential 


of the Milky Way. Peters (1975) adopted a model that takes 


into accout gas flowing along elliptic streamlines caused by 
a bar-like potential and showed that the general trend of 
the neutral hydrogen (HI) distribution is explained by the 
model. Liszt Burton] ( |1980| developed a tilted bar model 
that accounts for the observed HI distributions at different 
latitudes. 
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Observationally, Blitz & Spergel (1991) identified a bar¬ 


like distribution of stars in the near-infrared distribution b 
see also 


Matsumoto et al. 


(1982 


Hayakawa et al. (1981) 


Based on the above consideration and finding, |Binney et al.| 


(1991) showed that a bar-like potential naturally reproduces 


the noncircular motion by comparing observed and calcu¬ 
lated I — V (Galactic longitude - line-of-sight velocity from 
the Local Standard of Rest; LSR, hereafter) diagrams. The 
streaming motion driven by the bar potential provides an 
only viable explanation on the parallelogram, which has 
been explicitly discussed in literatures (Blitz 1994 |Koda 


Wada|2002| |Rodriguez-Fernandez Coml^ 2008 [Baba, 

Saitoh Wada||2010|[Molinari et al.|2011| ). 

We however find that the observed parallelogram in the 
2.6mm CO emission shows asymmetry in the positive and 


^ Infrared observation of the Galactic center region was also in¬ 


dependently done by Okuda et al. (1977). 
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negative velocities and in /, including CO features only at 
positive longitudes, I =3 degrees (clump 2) and 5.5 degrees 
(e.g., Bania|p^77 Torii et ^|201Qb| ). We also see that the 
central molecular zone (CMZ, hereafter) shows vertical CO 
features up to vertical distance ~ 200 pc (Enokiya et al. 
Torii et al.|2014a|b ), we discuss further later in 


2014 


this section. We note that these outstanding observed char¬ 
acteristics remain unexplored by the bar potential model. 

We note that few efforts have actually been made to 
elaborate detailed observed features of the parallelogram 
until recently by numerical simulations in the bar poten¬ 


tial. Rodriguez-Fernandez & Combes (2008) presented the 


first numerical simulations for the parallelogram. Their Fig¬ 
ure 14 shows a simulated parallelogram to be compared to 
the CMZ, while their result is yet crude at best with no 
detailed htting to the gas distribution in the l-v diagram. 
These authors concluded that the observed asymmetry of 
the CMZ cannot be due to lopsidedness of the bar and sug¬ 
gested that the asymmetry may be of an ad-hoc origin like 
infalling gas into the CMZ at /=1.3°. Their conclusion in¬ 
dicates that the bar potential alone cannot reproduce the 
asymmetry of the CMZ or the high-z CO features includ¬ 
ing the 1.3° complex. The difficulty in explaining the CMZ 
makes a contrast with the large-scale simulations in the bar 
potential over several kpc, which successfully reproduce the 
3-kpc arm and other major features outside the central kpc 
(e.g., Rodriguez-Fernandez Combes||2Q08 ). 

On the other hand, magnetic field is also supposed to 
play an important role in the Galactic bulge. Complicated 
structures, e.g., nonthermal filaments including Radio Arc, 
are observed in the Galactic center region, which probably 
attribute to magnetic fie lds (|Yusef-Zadeh, Morris Chaii^ 


1984 Tsuboi et al.|1986l|Chuss et al.|2003||Nishiyama et al. 


2010 Morris||2014| . |Crocker et al.| ( |201(5f ~gave a lower limit 

on the held strength, > 50 /xG, over the central 400 pc region 
from a non-thermal radio spectrum. By other considerations, 
it is argued that inferred magnetic held strength there is as 
strong as ~ mC ( Morris et al.|19^ Ferriere|2009 ). 

If such strong magnetic held is distributed in the bulge. 


it affects the global gas dynamics there (e.g., Sofue 2007). 
In fact, Fukui et al. (2006) discovered two large molecular 


loops, loop 1 and loop 2, at a radius of 700-pc from the 
center, which are triggered by magnetic buoyancy (Parker 
instability; Parker 1966); while the main focus in Parker 


(1966) was aimed at the Galactic disk, he already suggested 


that the magnetic activity might become even more im¬ 
portant in the central few 100 pc. Subsequently, |Fujishita| 
et al. (2009); Torii et al. ( 2010a|b ); Riquelme et al. ( |2010 ); 
Kudo et al. (2011) presented analyses of more detailed ob¬ 


servational properties of these molecular loops, and it was 
shown that the molecular loops are also well reproduced 
by magneto-hydrodynamical (MHD, hereafter) simulations 


(Shibata & Matsumoto 1991|[Machida et al.|2009|[Takahashi 


et al.||2009| ).' 

In spite of the broad recognition of the strong magnetic 
held in the Galactic center, there is little theoretical work 
that explores the role of the magnetic held in the dynamics 
of the gas component, except a limited number of attempts 


(e.g., Machida et al.|2QQ9 Kim Stone|2Ql'2 Machida et al. 


2013). We naturally expect that the magnetic held may play 


an essential role in the noncircular motion of the gas in the 



0.01 0.1 1 10 
R(kpc) 


Figure 1. Radial prohle of the equilibrium initial condition at 
the Galactic plane, top: Sound speed (dashed) and initial rota¬ 
tion speed (solid) at the Galactic plane. The rotation speed of 
the stellar component, Urot,* = (dotted) is also shown 

for comparison. Bottom: Initial density (dashed; left axis), gas 
pressure (solid; right axis), and magnetic pressure, pe = 

(dotted; right axis). pB is multiplied by 1000 times to ht in the 
panel. 


Galactic center region, which is the main aim of the present 
paper. 

In this paper, we investigate the evolution and the role 
of magnetic held in the central region of the Milky Way by 
a three-dimensional (3D) MHD simulation. We do not take 
into account a bar potential but an axisymmetric potential 
to focus on the role of the magnetic held in exciting the 
noncircular motion of gas. 


2 SIMULATION SETUP 


We treat the evolution of gas by solving MHD equations un¬ 
der an external axisymmetric Galactic gravitational poten¬ 
tial. We consider three components of gravitational sources: 
The supermassive blackhole (SMBH) at the Galactic cen¬ 
ter (component z = 1), a stellar bulge (i = 2), and a stel¬ 
lar disk {i = 3). Fo r the SMBH we assume a point mass . 
Ml = 4.4 X 10® Mq ( Genzel, Eisenhauer Gillessen||2010 ), 
where Mq is the solar mass. For the z = 2 and 3 components, 
we adopt a gravitational potential introduced by [Miyamotoj 


& Nagai (1975) for the bulge and disk components. The 


gravitation potential that includes these three components 
is written as 




-GMi 




( 1 ) 


where R and z are cylindrical radius and vertical distance, 
respectively, and the adopted values for Mi^ and bi are 
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Mi(lOiOM0) 

ai(kpc) 

^z(kpc) 

SMBH 1 

4.4 X 10-4 

0 

0 

Bulge 2 

2.05 

0 

0.495 

Disk 3 

25.47 

7.258 

0.52 


Table 1. Parameters for the gravitational potential. The param¬ 
eters for the super-massive black hole (SMBH) are from [Genzel,| 
[Eisenhauer &: Gillessen| ( |2010| ), and the parameters for the bulge 
and disk components are from [Miyamoto fc Nagaij ( |1975| >. 


summarized in Tableai = hi = 0 since we assume a point 
mass for the SMBH, and for the bulge and disk components 
we use the same values for Mi , ai , and hi in [Miyamoto fcj 
Nagai ( 1975| ) (see also Machida et ah (2009)). 

Assuming an equation of state for ideal gas, the gas 
pressure, 

p = pcl, ( 2 ) 

where p is gas density and Cs(oc a/T) is sound speed. We 
consider the spatial dependence of temperature (oc Cg) but 
assume that it is kept constant with time in each grid cell; 
the gas is assumed to be locally isothermal. Since we treat 
the Galactic bulge and disk by global simulation, the “sound 
speed” here more or less reflects the velocity dispersion of gas 
clouds rather than the actual temperature of each cloud. In 
the disk region, we assume Cg is comparable to the observed 
Bovy et al.|2012 ), 

(3) 


velocity dispersion (e.g. 


Cs,disk = 30km s , 
whereas this value is much smaller than the rotation speed 


ity dispersion, which is probably because of the noncircular 
motion, is obtained (e.g., Kent[p'9^ ). To mimic this larger 
velocity dispersion, we adopt 


Cs,bulge — O.GUrot,* — ^'^\/^ ’ 




(4) 


r(kpc) 

e 

d Nr Ng 

(0.01,60) 

(tt/O, 57r/6) 

(-7r,7r) 384 128 256 


Table 2. Simulation Box & Numerical Resolution 


Figure presents the radial distribution of the initial 
equilibrium profile at the Galactic plane. The rotation speed 
of the gas component shows a bump near R = 0.5 kpc, as¬ 
sociated with a peak of the sound speed there (top panel). 
In this region, because of the high temperature, the out¬ 
ward pressure-gradient force is not negligible, compared to 
the centrifugal force; the inward gravity is balanced by both 
the pressure-gradient force and the centrifugal force. To sat¬ 
isfy this radial force balance, the rotation speed should be 
kept considerably smaller than the rotation speed of the stel¬ 
lar component there to suppress the centrifugal force. This 
equilibrium configuration plays an important role in the am¬ 
plification of the magnetic field as will be discussed later. It 
should also be noted that the mass distribution of stars in 
the inner ~ 0.1 kpc has recently been obtained ( Sofue|2013 


Feldmeier et al.|2014| , which shows the existence of a certain 


amount of mass (~ 10®M©) there. We need to incorporate 
this contribution in future studies, which would give mod¬ 
erately faster rotation speed in the inner bulge. 

The bottom panel of Figure [^presents the initial den¬ 
sity, the gas pressure, and the magnetic pressure. The den¬ 
sity is expressed in units of number density, n cm“® = 
plpmu, with mean molecular weight, p = 1.2, where mu 
is the mass of a hydrogen atom. In this work we adopt the 
one-fluid approximation and do not distinguish ions, neu¬ 
tral atoms, and molecules, n with /i = 1.2 approximately 
corresponds to number density in units of hydrogen atoms. 
For rough estimates of molecular number density, which is 
dominated by H 2 , in dense clouds, we can use tt-Hs ~ 0.5n. 
Initially we set up weak vertical magnetic field with. 


B, = 0.71/iG 


R 


1 kpc 


( 8 ) 


where Urot,* is the rotation speed of the stellar component. 

We smoothly connect these two components to fix the 
radial dependence: 


2 


G,bulge 


1 — tanh 


fR-Rh\ 
V Ai^b ) 


+c; 


's,disk 


1 + tanh 


V Ai?b ). 


where we assume the boundary between the bulge and the 
disk, Rh = 1 kpc, with a smoothing width, Ai^b = 0.8 kpc. 

We start our simulation from an equilibrium configu¬ 
ration, in which the gravity is balanced with the pressure- 
gradient and centrifugal forces: 


in the region of > 0.035 kpc. A reason why we start from 
the weak magnetic field is that we would like to avoid the 
direct effect of the initial magnetic field on the dynamics of 
the gas. The initial magnetic pressure, pu = jSir, is less 
: (3) than < 10“® times the gas pressure as shown in the bottom 
panel of Figure namely the plasma /3 = pjpu > 10®. 

We update p, u, and B with time by solving ideal MHD 
equations, 

^+V-ipv) = 0, (9) 


1 dp 

pM 




and 


1 dp 
p dz 


dz 


= 0 . 


dv 


0, 

(6) 

^ dt 



and 


(7) 

II 


V)t,-V (p + g) + (£ • v) (10) 


= V X (v X B), 


( 11 ) 


We assume the initial gas density at the Galactic plane is 
proportional to the stellar density that fixes the gravita¬ 
tional potential, namely p oc Then, we determine the 

distribution of the initial density and rotation frequency, Q, 
to satisfy Equations ^ and 0- 


under the fixed gravitational potential (Equation [^. The 
numerical scheme we adopt is the second-order Godunov- 
GMoGGT method ( Sano, Inutsuka fc Miyama] [1999| ), in 
which we solve nonlinear Riemann problems with magnetic 
pressure at cell boundaries for compressive waves and adopt 
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Figure 2. Time evolution of the magnetic field strength averaged 
over the azimuthal and vertical directions (see text) at different 
radial locations, R = 0.1 kpc (solid), 0.5 kpc (dashed), and 2.0 
kpc (dotted). 


the consistent method of characteristics (CMoC) for the 
evolution of magnetic fields ( Clarke|[l996| [Stone Norman 


1992) under the constrained transport (CT) scheme ( [Evans 


& Hawley 1988). Since we assume the locally isothermal 
equation of state (Equation]^, we do not solve an energy 
equation. We perform our simulation in spherical coordi¬ 
nates^ (r, 0 ), instead of cylindrical coordinate, (i?, 0 , z), 

to resolve the central region by fine-scale grids. The size 
of the simulation box and the number of the grid points 
are summarized in Table In the 0 direction, the simula¬ 
tion domain covers di 7 r /3 from the Galactic plane; it does 
not cover the regions near the poles, and we need to set up 
an appropriate condition for the 0 boundaries, which is de¬ 
scribed below. The size of the radial grid, Ar, is enlarged in 
proportion to r. We analyze the numerical data mainly in 
the cylindrical coordinates by converting the primitive data 
in the spherical coordinates. 

We adopt the same boundary condition as in [Suzuki[ 


& Inutsuka (2014). We prescribe the outgoing condition for 


both mass and waves at the 0 boundaries, which was orig¬ 


inally developed for simulations for the solar wind ( Suzuki 
&: Inutsuka|2QQ5 2006), and the accretion condition at both 


the outer and inner r boundaries. Because of the outgoing 
boundaries at the 0 surfaces, the mass rapidly streams out 
by vertical outflows driven by effective magneto-turbulent 
pressure ( Suzuki Inutsuka|2009 2014). Therefore, the to¬ 
tal mass in the simulation box does not conserve but de¬ 
creases with time. This process probably overestimates the 
mass loss especially in the bulge region because our simula¬ 
tion box does not cover up to the sufficiently higher altitude 
( Suzuki, Muto Inutsuka|2010 Fromang et al.|20l3 ). In or- 



Figure 3. Snapshot views of the bulge region (R < 1 kpc) of 
the simulation ai t = 439.02 Myr (upper panel). Colours indi¬ 
cate density in units of number density, n cm“^, in logarithmic 
scale, white lines denote magnetic field lines. In the lower panel 
the central region is zoomed in and seen from a nearly edge-on 
angle to show ^-shaped field lines which are typical for Parker 
instability. The movies are available (see text and footnote). 


der to compensate the rapid mass decrease in the bulge, we 
use very large “density floor”, pfl of 0.8 times the local ini¬ 
tial value, pinit; if density becomes less than the floor value, 
Pfl = 0.8pinit, at each location, we recover it to pa. The den¬ 
sity floor is utilized in the only regions near the 0 surfaces 
where the mass is lost by vertical outflows. The density there 
is lower than that near the Galactic plane, and the effect of 
mass gain by this treatment is small and it only partially 
compensates the mass lost by vertical outflows. 

When we carry out the simulation, all the physical 
variables are non-dimensionalized by unit mass, munit = 

1 O^°M 0 , unit length, Gnit = 1 kpc, and unit time, tunit = 

. - 1/2 

= 5.15 Myr. We perform the simulation up 


/ Gruumt] 

\ unit / 


to t = SOyrtunit = 485.21 Myr. 
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X(kpc) 


Figure 4. Face-on views of density in units of n cm ^ (colour) and velocity field (arrows) at the Galactic plane at different times, 
t = 439.02 Myr (left) and 441.61 Myr (right). 



Figure 5. l-v diagrams at different times, t = 439.02 Myr (left) and 441.61 Myr (right). The grid points with high-density regions, 
n > 100 cm“^, in | 2 ;| < 0.2 kpc are used to derive the column density in units of 10^^cm“^ (colour) by integrating n along the direction 
of line of sight. In the horizontal axis, both Galactic longitude, I degree, (bottom axis) and the corresponding transverse distance, x kpc, 
at y = 0 kpc (top axis) are shown, where the solar system is located at (x, y) = (0, —8) kpc. 



I (deg) I (deg) 

Figure 6. Same as Figure but for a large I range. Regions with smaller density, n > 1 cm“^, are considered to derive the column 
density than for Figure 
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3 RESULTS 


3.1 Overview 


After the simulation starts, the amplification of the weak 
vertical seed magnetic held is initially triggered by magneto- 
rotational instability (MRI, hereafter; Velikhov|1959 Chan¬ 
drasekhar 1961 Balbus & Hawley 1991) and the vertical 
shear of the initial rotation prohle ( Suzuki Inutsuka|2014 


McNally Pessah|2014 ) to generate the radial and toroidal 


components. The wavelength of the most unstable mode of 
the MRI, which is captured by ~ one grid scale, is slightly 
underresolved at hrst because the initial held strength is 
weak. However, the wavelength of the most unstable mode, 
which is proportional to the held strength, becomes longer 
with the amplihcation of magnetic held, and it can be well 
resolved in the saturated state.The toroidal component is 
further amplihed by the winding due to the radial diher- 


ential rotation. Furthermore, magnetic buoyancy (Parker 
1966) also plays a role in amplifying the vertical compo¬ 


nent ( Nishikori, Machida Matsumoto|20Q6 Machida et al. 


2013). 


Figure presents the time evolution of the magnetic 
held strength averaged over azimuthal (0) and vertical [z) 
directions. 


\/W)Z = 




2n{z2 - Zi) 


( 12 ) 


at diherent radial locations, R = 0.1 kpc (solid), 0.5 kpc 
(dashed), and 2.0 kpc (dotted), where the vertical integra¬ 
tion is taken from zi — — 1 to Z 2 = 1 kpc. The growth 
of magnetic held is faster at smaller R because the growth 
rate of the MRI and the winding is scaled by rotation fre¬ 
quency, Q (see f 3.3), which decreases with R. Figure|^shows 
that the held strength saturates after t ~ 400 Myr at 

R = 0.1 and 0.5 kpc in the bulge region, while it is still 
gradually growing at R = 2 kpc until the end of the simula¬ 
tion, t — 307rtunit = 485.21 Myr. 

Figure [^presents 3D snapshots |^of the bulge region, 
R < 1 kpc, after the MHD turbulence is well developed 
at t = 439.02 Myr. The hgure exhibits turbulent gas and 
complicatedly tangled magnetic held lines in the bulge re¬ 
gion. One can see turbulent poloidal held lines excited by 
MRI and Parker instability, whereas the toroidal component 
slightly dominates the poloidal (radial and vertical) compo¬ 
nent as will be discussed in §3.3| In the lower panel that 
zooms in the central region, we pick up ^-shaped held lines 
excited by Parker instability, which are actually observed in 
the Galactic center region ( Fukui et al.||2006 ). This turbu¬ 
lent state is sustained in a quasi-steady manner after the 
held strength is saturated. In Figure we only show the 
evolution of the average held strength, and then one can¬ 
not see evolutionary properties of the magnetic conhgura- 
tion and turbulence. The movie of Figure]^ (see footnote for 
the link) demonstrates that the MRI and Parker instabil¬ 
ity continuously drive non-axysmmetric complex magnetic 
conhguration rather than developing coherent axisymmetric 
structure. Therefore, the system is in a quasi-steady satu- 


^ Movie of the simulation is available at 

WWW.ta.phys.nagoya-u.ac.jp/stakeru/research/galaxypot/glbdsk28inv_lagp2.mp4 


rated state in terms of magnetic turbulence, after the satu¬ 
ration of the held strength. 

Figure shows the velocity held with the density, n, at 
the Galactic plane {0 = 7r/2) at diherent times the left 
panel presents the result at the same time as in Figureand 
the right panel shows the result at 2.6 Myr later than the 
left panel, where this time diherence (= 2.6 Myr) roughly 
corresponds to 1/3 of the rotation time at R = 0.1 kpc. 
These two panels show that radial motion is stochastically 
excited particularly around R 0.5 kpc. Outward motion 
{vr > 0) appears to dominate, whereas one can also recog¬ 
nize inhows (vr < 0) in some regions. These behaviors are 
a result of the magnetic held as will be discussed in more 
detail later ( p.3| . 


3.2 l—v diagrams 

Figure]^ shows Galactic longitude-velocity (I — v) diagrams 
observed from the LSR at diherent times that correspond 
to Figure H The solar system is assumed to rot ate with 
240 km s“ in the clockwise direction at R = 8 kpc (Honma 
et al.|201^ and 0 = —'7r/2, namely {vx,Vy,Vz) = (—240, 0, 0) 


km s“^ at {x, y, z) = (0, —8, 0) kpc in the Gartesian coordi¬ 
nates. The contours indicate the column density in units of 
10^^cm“^, integrated along the direction of line of sight. We 
use the grid points with high-density regions where n > 100 
cm“^, because the high-density regions more or less trace 
the cool molecular gas, although our simulation treats one- 
huid gas and does not distinguish molecules, neutral atoms, 
or ions. 

Parallelogram-like shapes are observed in the central 
region as a consequence of the excited radial motion shown 
in Figure [4| G omparing the two panels, the shape changes 
with time |j on account of the time-dependent nature of 
the radial flows. Moreover, one can see that the shape is 
not symmetric because the distribution of the excited ra¬ 
dial flows is not axisymmetric or bisymmetric. Interestingly, 
the parallelogram shape obtained in the Milky way is not 
symmetric ( Dame, Hartmann Thaddeus|2001| Torii et al. 
|2010a[ then, the shape is not a “parallelogram” in a strict 
sense), which cannot be explained by a simple bar potential. 

Figure [^presents l-v diagrams for a larger |/| < 15 de¬ 
gree. We use regions with smaller density, n > 1 cm“^, to 
calculate surface density than for Figure in order to take 
into account the outer region where the density is smaller. 
These diagrams also show non-symmetric structure in the 
l-v plane, which indicates that non-axisymmetric gas distri¬ 
butions are extended to larger R. 

We examine the time-dependency of the parallelogram¬ 
like feature in Figure The solid lines in the top panel 
shows the time-evolution of the maximum and minimum 
velocities (r’max & ^’min) for column density > 0.5x lO^^cm”^ 
at / = 0 degree; they correspond to the positions where the 
upper and lower sides of the parallelogram cross the vertical 
line of / = 0 degree in the I — v diagram. The hgure shows 
that I’max and I’min vary within 50-150 km s“^ with time. 


^ Movie is available at 

WWW.ta.phys.nagoya-u.ac.jp/stakeru/research/galaxypot/faceon28inv_4.gif 

Movie is available at 

www.ta.phys.nagoya-u.ac.jp/stcikeru/research/galaxypot/pv-anim28inv_hr_2.gif 
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tinne(Myr) 



time(Myr) 


Figure 7. top: Time evolution of heliocentric velocities at Galac¬ 
tic longitude, I = 0 degree The solid lines indicate the max¬ 
imum and minimum velocities that give column density, E > 
0.5 X 10^^cm“^. The dashed lines trace the radial motion of the 
front-side (0 = —7r/2) and back-side {(p = tt/ 2) ai R = 0.5 kpc 
{{x,y) = (0, ^0.5) kpc in the Cartesian coordinates), bottom: 
Time evolution of the slope of the “parallelogram” near I = 0 
degree. The solid and dashed lines indicate the slope of the top 
and bottom sides, respectively. The slope is derived from the av¬ 
erage between I = ±0.72 degrees, which corresponds to ±100 
pc around the Galactic center. The shaded regions indicate the 
observational values, the region centered at 20 km s“^deg“^ for 
the top side of the parallelogram and the region centered at 55 
km s“^deg“^ for the bottom side. 


Basically I’max and I’min change independently of each other, 
and then, the shape of the symmetric parallelogram is not 
always kept but distorted with time as seen in Figure 

In Figure the time evolution of the motion of the gas 
at (x, y) = (0, ±0.5) kpc is also shown (dashed lines) for 
comparison. These lines indicate that the direction of the 
radial motion tends to be outward (r’ < 0 for the front side 
and i’ > 0 for the back side), whereas inward motions are 
also seen occasionally. 

© 2015 RAS, MNRAS 000 , [^?? 



Figure 8. Radial distribution of the plasma /3 (top) and the a 
value {bottom; solid line) which is the sum of the Maxwell (dashed 
line) and Reynolds (dotted line) stresses. The data are averaged 
over the azimuthal direction in full rotation (27r), over the vertical 
direction within | 2 ;| < 1 kpc, and the time t = 436.68-485.21 Myr. 


The bottom panel presents the slope of the parallelo¬ 
gram measured near / = 0, which is derived from the aver¬ 
age in the region with —0.72 < I < 0.72 deg, corresponding 
to X ~ ±100 pc. The solid and dashed lines denote the top 
and bottom sides of the parallelogram, which mostly reflect 
far and near sides with respect to the Galactic center, re¬ 
spectively. The top and bottom sides change independently 
because of the non-axisymmetric excitation of radial flows. 
The slopes vary from 15 to 90 km s“^deg“^, which covers 
the observed slope, ~ 20 km s“^deg“^, for the top side of 
the parallelogram and 55 km s“^deg“^ for the bottom 
side ( Torii et al.|2010a ). 


3.3 Origin of Noncircular Motion 
—Roles of Magnetic Field— 

Our MHD simulation shows that noncircular motions are 
excited, even though the axisymmetric gravitational poten¬ 
tial is adopted. Here we discuss the roles of the magnetic 
field in driving the radial gas flows by inspecting the numer¬ 
ical data. In Figures and we present the radial profile of 
various quantities averaged over azimuthal (0) and vertical 
(z) directions and over time t. We take the time-average of 
the final stage of the simulation from ti = 277rtunit = 436.68 
Myr to t 2 = 307rtunit = 485.21 Myr. Although the magnetic 
field is already amplified to achieve the saturated state in the 
bulge region, it is still developing outside the bulge (Figure 
1^. Therefore, in these figures we focus on the inner region. 
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R(kpc) 

Figure 9. Radial distribution of various quantities. The data are 
averaged over the azimuthal direction in full rotation (27r), over 
the vertical direction within | 2 ;| < 1 kpc (except for panel (a)), 
and the time t = 436.68-485.21 Myr. (see text), (a): The strength 
of differential rotation, d\nQ/d\nR at the Galactic plane (solid 
line); The average is not taken over the 2 : component but over 
the (j) component and time. The dashed line represents the initial 
condition. Rigid-body rotation {d\nQ,/d\n R = 0), ffat rotation 
(= —1), and Keplerian rotation (= —3/2) are shown by dotted 
lines for reference, (b): R (solid), ^(dashed), and 2 ; (dotted) com¬ 
ponents of root-mean squared (rms) magnetic field in units of /rG. 
The dotted line with a constant slope denotes the initial vertical 
magnetic field, (c): Root mean squared radial velocity, 

(solid), and mean radial velocity, (vr) (dashed and dotted lines 
for positive and negative values). 


R ^ 2 kpc. We take the average of a physical quantity, A, 
as follows: 

rt2 nZ2 

{A) = / dt d(j) dzA/[27T{t2 — ti){z2 — zi)], (13) 

J J — TT J Z\ 

where we set zi = — 1 and 2:2 = +1 kpc. For the quantities 
concerning velocity, we adopt density-weighted averages; the 
average of flow velocity and root-mean squared (rms) veloc¬ 
ity is taken by 

(.) . (.4) 

and 



Figure 10. Simplified schematic picture of the mechanism that 
drives radial fiows. Let us suppose three rings, ‘F, ‘ 2 ’, and ‘3’, that 
rotate with rotation frequency, lli, 112 , and Us, at the outer edge 
of each ring. If Hi > H 2 = H 3 , differential rotation is more in¬ 
tense in the red region (Region 2). The magnetic field is amplified 
more effectively there by the MRI and the field-line stretching, 
and the Maxwell stress, om, (Equation EH is larger than those 
in the neighboring zones, which gives the faster outward trans¬ 
port of the angular momentum in the red region. As a result. Hi 
decreases and H 2 increases. From the radial force balance, the in¬ 
ner edge of the red region moves inward because of the decrease 
of the centrifugal force, while the outer edge moves outward by 
the increase of the centrifugal force. Although the inhomogene¬ 
ity of angular momentum transport is distributed in a far more 
complicated fashion in our simulation, the essential point can be 
explained by this simple picture. 



while the simple average is taken for the quantities on mag¬ 
netic field. 

Figure [^presents the plasma /3 value. 
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Pb.3* 

Pb,2 ^ 1 r ^B, 3 

Figure 11. Additional mechanism that drives radial flows. The 
initial setting is the same as in Figure Radial magnetic held, 
which is triggered by MRI in our simulation, is amplifled to gen¬ 
eral toroidal held by the differential rotation selectively in the red 
region (Region 2). The magnetic pressure, pB, of fho toroidal held 
there is larger than pB in fho neighboring regions. Then, radial 
ffows are generated by the magnetic pressure gradient force. 


/«\ = AL = 

(PB> 


(16) 


which is the ratio of gas pressure to magnetic pressure (top 


panel), and the a value (bottom panel; solid line) ( [Shakura 
Sunyaev||1973| , 

(pvRSVff)) {BRB(f)) 


(a) = (qr) + (<aM) = 


{p) 


47r(p) 


(17) 


which is the sum of Reynolds (qr; dotted line) and Maxwell 
(qm; dashed line) stresses. The fluctuation component, 
is derived by subtracting the mean rotation. 


Svcj, = Vcf, - (v^). (18) 

The top panel of Figure shows that the plasma [3 is 
below 100 in the bulge region, < 1 kpc, which is much 
smaller than the initial value, 13 ~ 10^ — 10^, as a con¬ 
sequence of the amplification of the initial weak magnetic 
field. The bottom panel shows that a large a value > 10“^ 
is obtained. In the inner bulge, < 0.3 kpc, the Maxwell 
stress dominates the Reynolds stress. In the outer bulge, the 
Reynolds stress exceeds the Maxwell stress, and they both 
show bumpy structures, which reflects the complicated ro¬ 
tation profile as discussed from now. 

Figure]^ (a) presents the time and 0-averaged strength 
of the radial differential rotation, ^ In In R, at the Galac¬ 
tic plane. Note that d\nQ/d\nR = 0, —1, —3/2 correspond 


to rigid-body, fiat, and Keplerian rotations, respectively. The 
MRI is active when dhiQ/d\nR < 0, and its maximum 
growth rate to generate radial magnetic field is given by 


<^MRI,max _ 1 

Q “ 2 


d\nQ 


d\nR 


(19) 


( Balbus fc Hawley|1991 ). The amplification of toroidal field 
by winding in the radial differential rotation follows (e.g., 
Balbus Hawle^|1998 ) 


1 dB(j) _ ^ In Q 

Q dt ^ d\nR 


( 20 ) 


Equations (19) and (20) show that stronger differential ro¬ 


tation amplifies magnetic field faster by both MRI and field¬ 
line stretching. 

Figure [^a) shows that the initial rotation profile is 
more or less kept at the later time. The local minimum of 
d\nQ/d\nR is located at R ~ 0.6 kpc, where the rotation 
rapidly decreases with R because of the contribution from 
the gas pressure-gradient force (^. In this region, the mag¬ 
netic field is effectively amplified on account of the strong 
differential rotation. The magnetic field strength is ampli¬ 
fied to \B\ = -\- B‘^-\- Bz 0.1 — 1 mG in the bulge 

region, which is consistent with an observational lower limit, 
> 50 pG ( Grocker et al.|2010 ) and an empirical estimate, ~ 
mG, based on multiple observations ( Ferriere||2009 ). 

The effective amplification of the magnetic field around 
the location for the maximum differential rotation at R ~ 
0.6 kpc gives a bump in the Maxwell stress, am (Equation 
0 as shown by the dashed line in the bottom panel of 
Figure aM determines the outward transport of angular 
momentum (e.g., Lynden-Bell &: Pringle|1974 ) by magnetic 
field. The inhomogeneous distribution of am indicates that 
spatially dependent gain or loss of the angular momentum 
takes place, which is illustrated by a simplified cartoon in 
Figure in a region with the gain of angular momentum 
the rotation is accelerated, and vice versa. This leads directly 
to the increase or decrease of the centrifugal force, which 
excites radial motion. Although what is taking place in our 
simulation is more complicated, namely am is distributed in 
a non-axisymmetric manner, and both outward and inward 
flows are generated even at the same R as exhibited in Figure 
1^ this simple conceptual cartoon in Figure explains the 
essential point. 

In addition to the inhomogeneous transport of angular 
momentum, magnetic pressure also contributes to the exci¬ 
tation of radial flows. In the outer bulge region near R ^ 1 
kpc, the differential rotation is weak and the amplification 
of magnetic field is suppressed there. Figure |^b) exhibits 
the rapid decrease of the magnetic field strength with R in 
R ~ 0.5 kpc. The rapid decrease of B^ causes the radial 

pressure gradient, which drives outward radial 

flows (Figure pr]). 

We have discussed the two types of the processes, the 
inhomogeneous angular momentum transport and the mag¬ 
netic pressure-gradient force, that generate radial flows. By 


inspecting each term in the momentum equation (10) nu¬ 
merically, we found that the former dominates the latter by 
2 : 1 . 

Figure [^c) shows that the direction of the mean flow 
is outward (dashed line (vr) > 0) in 0.3 ~ R ~ 2 

kpc, which is expected from both processes. This radially 
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-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 

X (kpc) 

Figure 12. Density structure in the Galactic center region at 
t = 439.02 Myr. (Zoomed-in view of the left panel of Figure]^) 



0 0.2 0.4 0.6 0.8 1 

R (kpc) 


Figure 13. Vertical cut at0 = 7r/2ofl//3 (colours) and velocity 
field (arrow) at t = 439.02 Myr. 


outward flow also transports the angular momentum, which 
is a reason in part why the Reynolds stress (the bottom 
panel of Figure]^ shows a bump in 0.5 ~ R ~ 1 kpc. 
Compared to the mean flow, the rms velocity, (solid 

line) at R ~ 0.5 kpc is quite large, ~ 40 km s“^, because the 
MRI-triggered turbulence drives both inward and outward 
vr intermittently ( |Sano Inutsuka|2001 see also movies of 
Figures]^ and 1^, as discussed so far. This is the main reason 
why the simulated l-v diagrams show a thick parallelogram¬ 
like shape (Figure]^. 

3.4 Non-axisymmetric Structure 

As we have discussed in the previous subsection, the ve¬ 
locity field shows non-axisymmetric radial flows because of 
the intermittency of the MRI-triggered turbulence. The den¬ 
sity distribution also shows non-axisymmetric inhomogene¬ 


ity, which is typical for the MRI turbulence as well (Suzuki 
(V Inutsuka||2Q14 ), m the central region as exhibited in Fig¬ 
ure nr One can see multiple arm-like structure, and when 
measured at R = 50 pc, the density varies largely within, 
1500 < n < 3800 cm“^, which roughly corresponds to 
750 < nH 2 < 1900 cm“^ in molecular number density. 

It is well known that the Milky way possesses the cen¬ 
tral molecular zone (CMZ), which contains molecular gas 
with mass (5 — 10) x lO^M© and occupi es a volumetric 
filling factor ~ 0.1 within R < 200 pc ( Morris Ser- 
abyn||1996 ). The CMZ consists of non-axisymmetric clouds 
with complicated structure ( Oka et al.|2005 Kruijssen, Dale 
fc Longmore||2015 ). Our simulation implies that such non- 
axisymmetric clouds are a natural outcome of magnetic ac¬ 
tivity in the Galactic center. The molecular number density 
estimated from our simulation abov e is smaller than typica l 
observed values, ~ 10^ cm“^ (e.g., Nagayama et al.||2007 ). 
This is because our simulation assumes the one-fluid gas 
with the locally isothermal equation of state. In R < 200 
pc, the temperature is assumed to be constant ~ 10^ K 
in the simulation, which is much higher than the observed 
typical temperature, 30-200 K, of molecular clouds in the 


CMZ (e.g. Morris & Serabyn 1996). We need to treat the 
energy equation with heating and cooling to deal with dif¬ 
ferent phases of the gas material such as molecular clouds 
and warm neutral media that respectively possess different 
filling factors. The temperature of dense regions in the sim¬ 
ulation would be much lower if the radiative cooling was 
taken into account, and accordingly the density would be 
higher to satisfy the pressure balance. 


3.5 Outflows 

Figure presents the snapshot of a vertical cut at t = 
439.02 Myr, Here the colours indicate 1//5 = Pb/p; redder 
colours correspond to regions with relatively strong mag¬ 
netic pressure compared to gas pressure. The figure shows 
that flow patterns are quite complicated; one can see both 
inward flows to the bulge and outward streams from it. 
Although both inward and outward streams coexist, if we 
integrate the mass flux in space and time, the direction 
of the net gas flow is outward from the bulge; for exam¬ 
ple the mass loss rate measured at the spherical surface, 
r = 1.0 kpc, averaged during t = 436.68 — 485.21 Myr, gives 
M — j dt j d(j) j dO j drpVrr^sinO ^ 70 M© yr“^, and the 
kinetic energy luminosity, ~ ~ 10^® erg s It should 

be noted that this outflow rate might be suppressed when 
a stellar-bar potential was taken into account ( [Chandran, 
Cowley Morris|2000 ). 

In addition to radial flows, upward flows (winds) stream 
out of the 0 surfaces of the simulation box. The typical speed 
of these flows is several tens km s“^, and in some regions 
the speed exceeds 100 km s“^. Such vertical outflows possi¬ 
bly explain the observed vertical structure of molecular gas 
associated with the DHN ( Enokiya et al.|[2014 Torii et al. 
2014b). Recent observation further reveals CO emission at 


high Galactic latitude, 6, including several filamentary fea¬ 
tures in addition to diffuse extended halo-like CO gas up to 
2° in h above and below the CMZ (Torii et al. 2014a). A 


plausible interpretation is that these high-5 features are also 
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driven by the Parker instability; although they are filamen¬ 
tary without a clear loop-like shape having two foot points, 
the MHD numerical simulations by Machida et ah (2009) 
show that magnetic-hot at ion loops are generally not sym¬ 
metric and can be seen as an open single filament depending 
on the ambient density/field distribution and evolutionary 
effect. Further signature is associated with the double helix 
nebula (DHN, hereafter) toward (/,6) = (0.0 deg., 0.8 deg.) 
( Morris, Uchida Do|[20Q6 ), where a column of molecular 
gas of 200-pc height is found to be associated with the DHN 


which is nearly vertical to the plane (Enokiya et al. 2014 
Torii et al.|2014b ). Such a feature is explained as created by 


a magnetic tower formed above and below the central black 


hole according to preceding numerical simulations ( |Kato, 
Mineshige Shibata|2004 Machida et al.||2009 ). 


Recently, the existence of outflows or winds in the 


Galactic center region is extensively discussed ( Everett et al. 
2008 Crocker et al.|20lT ), and it might be closely linked to 


the formation mechanisms of the Eermi bubbles via mag¬ 
netic activity (Carretti et al. 2013). Our simulation is not 
directly applicable to the Eermi bubbles since the simula¬ 
tion does not cover the region near the polar axis. However, 
the above rough estimate of the energetics of the outflow 
in our simulation shows that magnetic activity is a reliable 
mechanism in driving outflows from the bulge. 


4 SUMMARY AND DISCUSSION 


We have investigated the excitation of radial gas flows in 
the Galactic center region by the 3D global MHD simula¬ 
tion with the axisymmetric gravitational potential. The ini¬ 
tial weak vertical magnetic field is amplihed by the MRI 
at the beginning, and in addition eventually by the Parker 
instability and the field-line stretching owing to the differ¬ 
ential rotation. In the bulge region, the final field strength 
is ~ 0.5 mG, and the plasma P ~ 1 — 40. Because of 

the gas pressure-gradient force, the equilibrium rotation fre¬ 
quency is non-monotonic with radial distance. The ampli¬ 
fication of the magnetic field is systematically more effec¬ 
tive in regions with stronger differential rotation. Then, the 
transport of the angular momentum is spatially dependent, 
and the acceleration or deceleration of the rotation speed 
occurs depending on regions. This breaks the radial pres¬ 
sure balance owing to the change of the centrifugal force, 
and excites radial flows of the gas in a time-dependent and 
non-axisymmetric manner. In addition, the radial compo¬ 
nent of the magnetic pressure gradient is produced, which 
also excites radial flows. As a result, the simulated l-v di¬ 
agram exhibits a time-dependent and asymmetric parallel¬ 
ogram shape. The rotation curve near the Galactic center 
exhibits complicated features ( Sofue|[2013 ). Our simulation 
implies that stochastically excited radial motion might be 
contaminated in these features. 

In order to interpret the parallelogram in the GMZ, the 
elliptical orbit in the stellar-bar potential has been used as 
the only viable model ( Binney et al.|19'^ ). The model may 
not be fully appropriate as the interpretation from an ob¬ 
servational point. Eirst, the bar potential alone is not able 
to create high-6 features up to 200 pc above the plane; these 
loops are better explained by the Parker instability ( Fukui 
et al.|2006 ). There is no other driving mechanism known to 


expel gas up to that height by only stellar gravity. Second, 
the parallelogram is not symmetric in the l-v diagram; the 
velocity gradient in the positive velocity, ~20 km s“^ deg“^, 
is significantly different by a factor of 3 from that in th e neg¬ 
ative velocity, ~55 km deg“^ ( Torii et 'aL]|2010a ). The 
present model shows that the velocity gradient of the paral¬ 
lelogram is naturally produced as a time dependent manner, 
whereas it is not clear how such a difference is explained in 
the bar potential model. 

On the other hand, we see some features that are not re¬ 
produced in the present simulation. The large broad features 
like the clump 2 and the 5.5 deg feature are not reproduced 
in the present simulation. In our simulation we have as¬ 
sumed one-fluid gas and neglected heating and cooling with 
a simplified treatment of the locally isothermal gas. In or¬ 
der to treat dense molecular gas in a quantitative fashion, 
we should replace these simplifications. If radiative cooling, 
as well as adiabatic heating and cooling, is included, large 
density contrast between cool molecular gas and warm/hot 
gas will be naturally reproduced, which may explain the ob¬ 
served fine-scale features. 

We have shown that the magnetic activity possibly 
produces the overall trend of the observed asymmetric 
parallelogram-shape in the I — v diagram even though the 
axisymmetric gravitational potential is considered. It is also 
desirable to test quantitatively the importance of the mag¬ 
netic activity in contrast to role of the stellar bar potential in 
the noncircular motion of the gas after the cooling/heating 
effect is taken into account in our simulation. 

Although in this paper we have focused on the bulge 
region of our Milky Way, results of our simulation can be 
applied to other galaxies. For example, [Sakamoto et al.| 
(2006) observed NGG 253 by the Submillimeter Array and 
reported that they found an expanding circumnuclear disk 
with ~ 50 km s“^. The obtained position velocity diagrams 
show asymmetric parallelogram features in different wave¬ 
lengths. These features are qualitatively similar to but quan¬ 
titatively different from those obtained in the Milky Way. A 
possible explanation of the difference is the time variability 
inhering in the MHD turbulence, based on our simulation. 
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